function fid = fid2syn(varargin)

% FIDSYN synthesizes FID from complex exponential time-domain model.  
% FIDSYN assumes that time samples are at [0:NN-1]. Time-scaling must be
% accomplished by scaling the input parameters appropriately.
%
% FIDSYN(MODEL,KBASIS,NN)
%
% MODEL describes the time-domain model at each voxel.  
% If MODEL is a structure:
%    MODEL.amp is an image array of complex magnitude image Mxy.
%    MODEL.exp is an image array of decay+j*frequency parameter exp(A).
% If MODEL is a vector:
%    MODEL(1:end/2) is complex magnitude vector Mxy.  
%    MODEL(end/2+1:end) is A in decay+j*frequency vector exp(A).
%
% If KBASIS is a structure:
%    KBASIS.EKXX is exp(kx(1:NN)'*(2j*pi*xx)), where kx are kx coordinates and
%    x are x coordinates
%    KBASIS.EKYY is exp(ky(1:NN)'*(2j*pi*yy)), where ky are ky coordinates and
%    y are y coordinates
%    KBASIS.IX is a vector of column coordinates of the FOV
%    KBASIS.IY is a vector of row coordinates of the FOV
%    KBASIS.IXX is a vector of indices to step through one dimension of
%    coordinates of the FOV
% 
% If KBASIS is a matrix, then KBASIS is the array of ideal local
% responses to the encoding gradients for unit-magnitude on-resonance
% magnetization vectors at each pixel (x,y) location.  These are "basis"
% functions of a sort:   kbasis=exp(-j*2*pi*(kx(1:NF)*x+ky(1:NF)*y));
%
% NN is the desired length of the calculated FID Y.
%
% FIDSYN(MODEL,KX,KY,NN) computes the FID directly from the model and the
% kx and ky coordinates.  It assumes that MODEL is in IM structure.

model = varargin{1};
if nargin==3
    kbasis = varargin{2};
    NN = varargin{3};
else
    kx = varargin{2};
    ky = varargin{3};
    NN = varargin{4};
    kbasis = 0;
end
wfr = model.wfr;
wfi = model.wfi;
fid = zeros(1,NN,'single');

if isstruct(model)
    if isfield(model,'N')
        NM = model.N;
    else
        NM = numel(model.amp);
    end
else
    NM = length(model)/2;  % NM is the number of voxels in the model
end
indv1 = single(NM+1:2*NM);
indv2 = single(2*NM+1:3*NM);

if isstruct(kbasis)
    % use separable calculation
    at = zeros(length(kbasis.ix1d));
    bt = zeros(length(kbasis.ix1d));
    vec_ind = kbasis.ix+model.reso*(kbasis.iy-1);
    if isstruct(model)
        model = im2vec_model2(model);
        at(vec_ind) = model.pvec(indv1);
        bt(vec_ind) = weight_cexp(model.pvec(indv2),wfr+1i*wfi);
    else
        at(vec_ind) = model(indv1);
        bt(vec_ind) = weight_cexp(model(indv2),wfr+1i*wfi);
    end
    bo = at;
    fid(1) = kbasis.ekxx(1,:)*bo*kbasis.ekyy(1,:).';
    for nr = 2:NN;
        bo = bo.*bt;
        fid(nr) = kbasis.ekxx(nr,:)*bo*kbasis.ekyy(nr,:).';
    end
else
    % extract parameters
    if isstruct(model)
        % vectorize images
        model = im2vec_model2(model);
        at = model.pvec(indv1);
        bt = weight_cexp(model.pvec(indv2),wfr+1i*wfi);
    else
        at = model(indv1);
        bt = weight_cexp(model(indv2),wfr+1i*wfi);
    end
    % calculate FID via exponential decay and frequency response basis
    bo = at;
    if nargin==3
        fid(1) = at*kbasis(:,1);
        for nr = 2:NN
            bo = bo.*bt;
            fid(nr) = bo*kbasis(:,nr);
        end
    else
        fid(1) = at*exp(-1i*2*pi*(kx(1)*model.x(:)+ky(1)*model.y(:)));
        for nr = 2:NN
            bo = bo.*bt;
            fid(nr) = bo*exp(-1i*2*pi*(kx(nr)*model.x(:)+ky(nr)*model.y(:)));
        end
    end
end